home *** CD-ROM | disk | FTP | other *** search
/ STraTOS 1997 April & May / STraTOS 1 - 1997 April & May.iso / CD01 / INTERNET / SITES / LITTLE / P3SRC.ZIP / ATARI / SUPER.C < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-15  |  31.3 KB  |  1,640 lines
  1. /****************************************************************************
  2. *                   super.c
  3. *
  4. *  This module implements functions that manipulate superellipsoids.
  5. *
  6. *  Original code written by Alexander Enzmann.
  7. *  Adaption to POV-Ray by Dieter Bayer [DB].
  8. *
  9. *  from Persistence of Vision(tm) Ray Tracer
  10. *  Copyright 1996 Persistence of Vision Team
  11. *---------------------------------------------------------------------------
  12. *  NOTICE: This source code file is provided so that users may experiment
  13. *  with enhancements to POV-Ray and to port the software to platforms other
  14. *  than those supported by the POV-Ray Team.  There are strict rules under
  15. *  which you are permitted to use this file.  The rules are in the file
  16. *  named POVLEGAL.DOC which should be distributed with this file. If
  17. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  18. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  19. *  Forum.  The latest version of POV-Ray may be found there as well.
  20. *
  21. * This program is based on the popular DKB raytracer version 2.12.
  22. * DKBTrace was originally written by David K. Buck.
  23. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  24. *
  25. *****************************************************************************/
  26.  
  27. /****************************************************************************
  28. *
  29. *  Explanation:
  30. *
  31. *    Superellipsoids are defined by the implicit equation
  32. *
  33. *      f(x,y,z) = (|x|^(2/e) + |y|^(2/e))^(e/n) + |z|^(2/n) - 1
  34. *
  35. *    Where e is the east/west exponent and n is the north/south exponent.
  36. *
  37. *    NOTE: The exponents are precomputed and stored in the Power element.
  38. *
  39. *    NOTE: Values of e and n that are close to degenerate (e.g.,
  40. *          below around 0.1) appear to give the root solver fits.
  41. *          Not sure quite where the problem lays just yet.
  42. *
  43. *  Syntax:
  44. *
  45. *    superellipsoid { <e, n> }
  46. *
  47. *
  48. *  ---
  49. *
  50. *  Oct 1994 : Creation.
  51. *
  52. *****************************************************************************/
  53.  
  54. #include "frame.h"
  55. #include "povray.h"
  56. #include "vector.h"
  57. #include "povproto.h"
  58. #include "bbox.h"
  59. #include "matrices.h"
  60. #include "objects.h"
  61. #include "super.h"
  62.  
  63. /*****************************************************************************
  64. * Local preprocessor defines
  65. ******************************************************************************/
  66.  
  67. /* Minimal intersection depth for a valid intersection. */
  68.  
  69. #define DEPTH_TOLERANCE 1.0e-4
  70.  
  71. /* If |x| < ZERO_TOLERANCE it is regarded to be 0. */
  72.  
  73. #define ZERO_TOLERANCE 1.0e-10
  74.  
  75. /* This is not the signum function because SGNX(0) is 1. */
  76.  
  77. #define SGNX(x) (((x) >= 0.0) ? 1.0 : -1.0)
  78.  
  79. #define MIN_VALUE -1.0
  80. #define MAX_VALUE  1.0
  81.  
  82. #define MAX_ITERATIONS 20
  83.  
  84. #define PLANECOUNT 9
  85.  
  86. /*****************************************************************************
  87. * Local typedefs
  88. ******************************************************************************/
  89.  
  90. /*****************************************************************************
  91. * Static functions
  92. ******************************************************************************/
  93.  
  94. static int intersect_superellipsoid PARAMS((RAY *Ray, SUPERELLIPSOID *Superellipsoid, ISTACK *Depth_Stack));
  95. static int intersect_box PARAMS((VECTOR P, VECTOR D, DBL *dmin, DBL *dmax));
  96. static DBL power PARAMS((DBL x, DBL e));
  97. static DBL evaluate_superellipsoid PARAMS((VECTOR P, SUPERELLIPSOID *Superellipsoid));
  98. static int compdists PARAMS((CONST void *in_a, CONST void *in_b));
  99.  
  100. static int find_ray_plane_points PARAMS((VECTOR P, VECTOR D, int cnt, DBL *dists, DBL mindist, DBL maxdist));
  101. static void solve_hit1 PARAMS((SUPERELLIPSOID *Super, DBL v0, VECTOR tP0, DBL v1, VECTOR tP1, VECTOR P));
  102. static int check_hit2 PARAMS((SUPERELLIPSOID *Super, VECTOR P, VECTOR D, DBL t0, VECTOR P0, DBL v0, DBL t1, DBL *t, VECTOR Q));
  103. static int insert_hit PARAMS((OBJECT *Object, RAY *Ray, DBL Depth, ISTACK *Depth_Stack));
  104.  
  105. static int   All_Superellipsoid_Intersections PARAMS((OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack));
  106. static int   Inside_Superellipsoid PARAMS((VECTOR point, OBJECT *Object));
  107. static void  Superellipsoid_Normal PARAMS((VECTOR Result, OBJECT *Object, INTERSECTION *Inter));
  108. static void  *Copy_Superellipsoid PARAMS((OBJECT *Object));
  109. static void  Translate_Superellipsoid PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  110. static void  Rotate_Superellipsoid PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  111. static void  Scale_Superellipsoid PARAMS((OBJECT *Object, VECTOR Vector, TRANSFORM *Trans));
  112. static void  Transform_Superellipsoid PARAMS((OBJECT *Object, TRANSFORM *Trans));
  113. static void  Invert_Superellipsoid PARAMS((OBJECT *Object));
  114. static void  Destroy_Superellipsoid PARAMS((OBJECT *Object));
  115.  
  116. /*****************************************************************************
  117. * Local variables
  118. ******************************************************************************/
  119.  
  120. static METHODS Superellipsoid_Methods =
  121. {
  122.   All_Superellipsoid_Intersections,
  123.   Inside_Superellipsoid, Superellipsoid_Normal,
  124.   Copy_Superellipsoid,
  125.   Translate_Superellipsoid, Rotate_Superellipsoid,
  126.   Scale_Superellipsoid, Transform_Superellipsoid, Invert_Superellipsoid, Destroy_Superellipsoid
  127. };
  128.  
  129. static DBL planes[PLANECOUNT][4] =
  130.   {{1, 1, 0, 0}, {1,-1, 0, 0},
  131.    {1, 0, 1, 0}, {1, 0,-1, 0},
  132.    {0, 1, 1, 0}, {0, 1,-1, 0},
  133.    {1, 0, 0, 0}, 
  134.    {0, 1, 0, 0}, 
  135.    {0, 0, 1, 0}};
  136.  
  137. /*****************************************************************************
  138. *
  139. * FUNCTION
  140. *
  141. *   All_Superellipsoid_Intersections
  142. *
  143. * INPUT
  144. *
  145. *   Object      - Object
  146. *   Ray         - Ray
  147. *   Depth_Stack - Intersection stack
  148. *   
  149. * OUTPUT
  150. *
  151. *   Depth_Stack
  152. *   
  153. * RETURNS
  154. *
  155. *   int - TRUE, if a intersection was found
  156. *   
  157. * AUTHOR
  158. *
  159. *   Dieter Bayer
  160. *   
  161. * DESCRIPTION
  162. *
  163. *   Determine ray/superellipsoid intersection and clip intersection found.
  164. *
  165. * CHANGES
  166. *
  167. *   Oct 1994 : Creation.
  168. *
  169. ******************************************************************************/
  170.  
  171. static int All_Superellipsoid_Intersections(Object, Ray, Depth_Stack)
  172. OBJECT *Object;
  173. RAY *Ray;
  174. ISTACK *Depth_Stack;
  175. {
  176.   Increase_Counter(stats[Ray_Superellipsoid_Tests]);
  177.  
  178.   if (intersect_superellipsoid(Ray, (SUPERELLIPSOID *)Object, Depth_Stack))
  179.   {
  180.     Increase_Counter(stats[Ray_Superellipsoid_Tests_Succeeded]);
  181.  
  182.     return(TRUE);
  183.   }
  184.   else
  185.   {
  186.     return(FALSE);
  187.   }
  188. }
  189.  
  190.  
  191.  
  192. /*****************************************************************************
  193. *
  194. * FUNCTION
  195. *
  196. *   intersect_superellipsoid
  197. *
  198. * INPUT
  199. *
  200. *   Ray            - Ray
  201. *   Superellipsoid - Superellipsoid
  202. *   Depth_Stack    - Depth stack
  203. *   
  204. * OUTPUT
  205. *
  206. *   Intersection
  207. *   
  208. * RETURNS
  209. *
  210. *   int - TRUE if intersections were found.
  211. *   
  212. * AUTHOR
  213. *
  214. *   Alexander Enzmann, Dieter Bayer
  215. *   
  216. * DESCRIPTION
  217. *
  218. *   Determine ray/superellipsoid intersection.
  219. *
  220. * CHANGES
  221. *
  222. *   Oct 1994 : Creation.
  223. *
  224. ******************************************************************************/
  225.  
  226. static int intersect_superellipsoid(Ray, Superellipsoid, Depth_Stack)
  227. RAY *Ray;
  228. SUPERELLIPSOID *Superellipsoid;
  229. ISTACK *Depth_Stack;
  230. {
  231.   int i, cnt, Found = FALSE;
  232.   DBL dists[PLANECOUNT+2];
  233.   DBL t, t1, t2, v0, v1, len;
  234.   VECTOR P, D, P0, P1, P2, P3;
  235.  
  236.   /* Transform the ray into the superellipsoid space. */
  237.  
  238.   MInvTransPoint(P, Ray->Initial, Superellipsoid->Trans);
  239.  
  240.   MInvTransDirection(D, Ray->Direction, Superellipsoid->Trans);
  241.  
  242.   VLength(len, D);
  243.  
  244.   VInverseScaleEq(D, len);
  245.  
  246.   /* Intersect superellipsoid's bounding box. */
  247.  
  248.   if (!intersect_box(P, D, &t1, &t2))
  249.   {
  250.     return(FALSE);
  251.   }
  252.  
  253.   /* Test if superellipsoid lies 'behind' the ray origin. */
  254.  
  255.   if (t2 < DEPTH_TOLERANCE)
  256.   {
  257.     return(FALSE);
  258.   }
  259.  
  260.   cnt = 0;
  261.  
  262.   if (t1 < DEPTH_TOLERANCE)
  263.   {
  264.     t1 = DEPTH_TOLERANCE;
  265.   }
  266.  
  267.   dists[cnt++] = t1;
  268.   dists[cnt++] = t2;
  269.  
  270.   /* Intersect ray with planes cutting superellipsoids in pieces. */
  271.  
  272.   cnt = find_ray_plane_points(P, D, cnt, dists, t1, t2);
  273.  
  274.   if (cnt <= 1)
  275.   {
  276.     return(FALSE);
  277.   }
  278.  
  279.   VEvaluateRay(P0, P, dists[0], D)
  280.  
  281.   v0 = evaluate_superellipsoid(P0, Superellipsoid);
  282.  
  283.   if (fabs(v0) < ZERO_TOLERANCE)
  284.   {
  285.     if (insert_hit((OBJECT *)Superellipsoid, Ray, dists[0] / len, Depth_Stack))
  286.     {
  287.       if (Superellipsoid->Type & IS_CHILD_OBJECT)
  288.       {
  289.         Found = TRUE;
  290.       }
  291.       else
  292.       {
  293.         return(TRUE);
  294.       }
  295.     }
  296.   }
  297.  
  298.   for (i = 1; i < cnt; i++)
  299.   {
  300.     VEvaluateRay(P1, P, dists[i], D)
  301.  
  302.     v1 = evaluate_superellipsoid(P1, Superellipsoid);
  303.  
  304.     if (fabs(v1) < ZERO_TOLERANCE)
  305.     {
  306.       if (insert_hit((OBJECT *)Superellipsoid, Ray, dists[i] / len, Depth_Stack))
  307.       {
  308.         if (Superellipsoid->Type & IS_CHILD_OBJECT)
  309.         {
  310.           Found = TRUE;
  311.         }
  312.         else
  313.         {
  314.           return(TRUE);
  315.         }
  316.       }
  317.     }
  318.     else
  319.     {
  320.       if (v0 * v1 < 0.0)
  321.       {
  322.         /* Opposite signs, there must be a root between */
  323.  
  324.         solve_hit1(Superellipsoid, v0, P0, v1, P1, P2);
  325.  
  326.         VSub(P3, P2, P);
  327.  
  328.         VLength(t, P3);
  329.  
  330.         if (insert_hit((OBJECT *)Superellipsoid, Ray, t / len, Depth_Stack))
  331.         {
  332.           if (Superellipsoid->Type & IS_CHILD_OBJECT)
  333.           {
  334.             Found = TRUE;
  335.           }
  336.           else
  337.           {
  338.             return(TRUE);
  339.           }
  340.         }
  341.       }
  342.       else
  343.       {
  344.         /* 
  345.          * Although there was no sign change, we may actually be approaching
  346.          * the surface. In this case, we are being fooled by the shape of the
  347.          * surface into thinking there isn't a root between sample points. 
  348.          */
  349.  
  350.         if (check_hit2(Superellipsoid, P, D, dists[i-1], P0, v0, dists[i], &t, P2))
  351.         {
  352.           if (insert_hit((OBJECT *)Superellipsoid, Ray, t / len, Depth_Stack))
  353.           {
  354.             if (Superellipsoid->Type & IS_CHILD_OBJECT)
  355.             {
  356.               Found = TRUE;
  357.             }
  358.             else
  359.             {
  360.               return(TRUE);
  361.             }
  362.           }
  363.           else
  364.           {
  365.             break;
  366.           }
  367.         }
  368.       }
  369.     }
  370.  
  371.     v0 = v1;
  372.  
  373.     Assign_Vector(P0, P1);
  374.   }
  375.  
  376.   return(Found);
  377. }
  378.  
  379.  
  380.  
  381. /*****************************************************************************
  382. *
  383. * FUNCTION
  384. *
  385. *   Inside_Superellipsoid
  386. *
  387. * INPUT
  388. *
  389. *   IPoint - Intersection point
  390. *   Object - Object
  391. *   
  392. * OUTPUT
  393. *   
  394. * RETURNS
  395. *
  396. *   int - TRUE if inside
  397. *   
  398. * AUTHOR
  399. *
  400. *   Dieter Bayer
  401. *   
  402. * DESCRIPTION
  403. *
  404. *   Test if a point lies inside the superellipsoid.
  405. *
  406. * CHANGES
  407. *
  408. *   Oct 1994 : Creation.
  409. *
  410. ******************************************************************************/
  411.  
  412. static int Inside_Superellipsoid(IPoint, Object)
  413. VECTOR IPoint;
  414. OBJECT *Object;
  415. {
  416.   DBL val;
  417.   VECTOR P;
  418.   SUPERELLIPSOID *Superellipsoid = (SUPERELLIPSOID *)Object;
  419.  
  420.   /* Transform the point into the superellipsoid space. */
  421.  
  422.   MInvTransPoint(P, IPoint, Superellipsoid->Trans);
  423.  
  424.   val = evaluate_superellipsoid(P, Superellipsoid);
  425.  
  426.   if (val < EPSILON)
  427.   {
  428.     return(!Test_Flag(Superellipsoid, INVERTED_FLAG));
  429.   }
  430.   else
  431.   {
  432.     return(Test_Flag(Superellipsoid, INVERTED_FLAG));
  433.   }
  434. }
  435.  
  436.  
  437.  
  438. /*****************************************************************************
  439. *
  440. * FUNCTION
  441. *
  442. *   Superellipsoid_Normal
  443. *
  444. * INPUT
  445. *
  446. *   Result - Normal vector
  447. *   Object - Object
  448. *   Inter  - Intersection found
  449. *   
  450. * OUTPUT
  451. *
  452. *   Result
  453. *   
  454. * RETURNS
  455. *   
  456. * AUTHOR
  457. *
  458. *   Dieter Bayer
  459. *   
  460. * DESCRIPTION
  461. *
  462. *   Calculate the normal of the superellipsoid in a given point.
  463. *
  464. * CHANGES
  465. *
  466. *   Oct 1994 : Creation.
  467. *
  468. ******************************************************************************/
  469.  
  470. static void Superellipsoid_Normal(Result, Object, Inter)
  471. OBJECT *Object;
  472. VECTOR Result;
  473. INTERSECTION *Inter;
  474. {
  475.   DBL k, x, y, z;
  476.   VECTOR P, N, E;
  477.   SUPERELLIPSOID *Superellipsoid = (SUPERELLIPSOID *)Object;
  478.  
  479.   /* Transform the point into the superellipsoid space. */
  480.  
  481.   MInvTransPoint(P, Inter->IPoint, Superellipsoid->Trans);
  482.  
  483.   Assign_Vector(E, Superellipsoid->Power);
  484.  
  485.   x = fabs(P[X]);
  486.   y = fabs(P[Y]);
  487.   z = fabs(P[Z]);
  488.  
  489.   k = power(power(x, E[X]) + power(y, E[X]), E[Y] - 1.0);
  490.  
  491.   N[X] = k * SGNX(P[X]) * power(x, E[X] - 1.0);
  492.   N[Y] = k * SGNX(P[Y]) * power(y, E[X] - 1.0);
  493.   N[Z] =     SGNX(P[Z]) * power(z, E[Z] - 1.0);
  494.  
  495.   /* Transform the normalt out of the superellipsoid space. */
  496.  
  497.   MTransNormal(Result, N, Superellipsoid->Trans);
  498.  
  499.   VNormalize(Result, Result);
  500. }
  501.  
  502.  
  503.  
  504. /*****************************************************************************
  505. *
  506. * FUNCTION
  507. *
  508. *   Translate_Superellipsoid
  509. *
  510. * INPUT
  511. *
  512. *   Object - Object
  513. *   Vector - Translation vector
  514. *
  515. * OUTPUT
  516. *
  517. *   Object
  518. *
  519. * RETURNS
  520. *
  521. * AUTHOR
  522. *
  523. *   Dieter Bayer
  524. *
  525. * DESCRIPTION
  526. *
  527. *   Translate a superellipsoid.
  528. *
  529. * CHANGES
  530. *
  531. *   Oct 1994 : Creation.
  532. *
  533. ******************************************************************************/
  534.  
  535. static void Translate_Superellipsoid(Object, Vector, Trans)
  536. OBJECT *Object;
  537. VECTOR Vector;
  538. TRANSFORM *Trans;
  539. {
  540.   Transform_Superellipsoid(Object, Trans);
  541. }
  542.  
  543.  
  544.  
  545. /*****************************************************************************
  546. *
  547. * FUNCTION
  548. *
  549. *   Rotate_Superellipsoid
  550. *
  551. * INPUT
  552. *
  553. *   Object - Object
  554. *   Vector - Rotation vector
  555. *   
  556. * OUTPUT
  557. *
  558. *   Object
  559. *   
  560. * RETURNS
  561. *   
  562. * AUTHOR
  563. *
  564. *   Dieter Bayer
  565. *   
  566. * DESCRIPTION
  567. *
  568. *   Rotate a superellipsoid.
  569. *
  570. * CHANGES
  571. *
  572. *   Oct 1994 : Creation.
  573. *
  574. ******************************************************************************/
  575.  
  576. static void Rotate_Superellipsoid(Object, Vector, Trans)
  577. OBJECT *Object;
  578. VECTOR Vector;
  579. TRANSFORM *Trans;
  580. {
  581.   Transform_Superellipsoid(Object, Trans);
  582. }
  583.  
  584.  
  585.  
  586. /*****************************************************************************
  587. *
  588. * FUNCTION
  589. *
  590. *   Scale_Superellipsoid
  591. *
  592. * INPUT
  593. *
  594. *   Object - Object
  595. *   Vector - Scaling vector
  596. *   
  597. * OUTPUT
  598. *
  599. *   Object
  600. *   
  601. * RETURNS
  602. *   
  603. * AUTHOR
  604. *
  605. *   Dieter Bayer
  606. *   
  607. * DESCRIPTION
  608. *
  609. *   Scale a superellipsoid.
  610. *
  611. * CHANGES
  612. *
  613. *   Oct 1994 : Creation.
  614. *
  615. ******************************************************************************/
  616.  
  617. static void Scale_Superellipsoid(Object, Vector, Trans)
  618. OBJECT *Object;
  619. VECTOR Vector;
  620. TRANSFORM *Trans;
  621. {
  622.   Transform_Superellipsoid(Object, Trans);
  623. }
  624.  
  625.  
  626.  
  627. /*****************************************************************************
  628. *
  629. * FUNCTION
  630. *
  631. *   Transform_Superellipsoid
  632. *
  633. * INPUT
  634. *
  635. *   Object - Object
  636. *   Trans  - Transformation to apply
  637. *   
  638. * OUTPUT
  639. *
  640. *   Object
  641. *   
  642. * RETURNS
  643. *   
  644. * AUTHOR
  645. *
  646. *   Dieter Bayer
  647. *   
  648. * DESCRIPTION
  649. *
  650. *   Transform a superellipsoid and recalculate its bounding box.
  651. *
  652. * CHANGES
  653. *
  654. *   Oct 1994 : Creation.
  655. *
  656. ******************************************************************************/
  657.  
  658. static void Transform_Superellipsoid(Object, Trans)
  659. OBJECT *Object;
  660. TRANSFORM *Trans;
  661. {
  662.   Compose_Transforms(((SUPERELLIPSOID *)Object)->Trans, Trans);
  663.  
  664.   Compute_Superellipsoid_BBox((SUPERELLIPSOID *)Object);
  665. }
  666.  
  667.  
  668.  
  669. /*****************************************************************************
  670. *
  671. * FUNCTION
  672. *
  673. *   Invert_Superellipsoid
  674. *
  675. * INPUT
  676. *
  677. *   Object - Object
  678. *   
  679. * OUTPUT
  680. *
  681. *   Object
  682. *   
  683. * RETURNS
  684. *   
  685. * AUTHOR
  686. *
  687. *   Dieter Bayer
  688. *   
  689. * DESCRIPTION
  690. *
  691. *   Invert a superellipsoid.
  692. *
  693. * CHANGES
  694. *
  695. *   Oct 1994 : Creation.
  696. *
  697. ******************************************************************************/
  698.  
  699. static void Invert_Superellipsoid(Object)
  700. OBJECT *Object;
  701. {
  702.   Invert_Flag(Object, INVERTED_FLAG);
  703. }
  704.  
  705.  
  706.  
  707. /*****************************************************************************
  708. *
  709. * FUNCTION
  710. *
  711. *   Create_Superellipsoid
  712. *
  713. * INPUT
  714. *   
  715. * OUTPUT
  716. *   
  717. * RETURNS
  718. *
  719. *   SUPERELLIPSOID * - new superellipsoid
  720. *   
  721. * AUTHOR
  722. *
  723. *   Dieter Bayer
  724. *   
  725. * DESCRIPTION
  726. *
  727. *   Create a new superellipsoid.
  728. *
  729. * CHANGES
  730. *
  731. *   Oct 1994 : Creation.
  732. *
  733. ******************************************************************************/
  734.  
  735. SUPERELLIPSOID *Create_Superellipsoid()
  736. {
  737.   SUPERELLIPSOID *New;
  738.  
  739.   New = (SUPERELLIPSOID *)POV_MALLOC(sizeof(SUPERELLIPSOID), "superellipsoid");
  740.  
  741.   INIT_OBJECT_FIELDS(New,SUPERELLIPSOID_OBJECT,&Superellipsoid_Methods)
  742.  
  743.   New->Trans = Create_Transform();
  744.  
  745.   Make_Vector(New->Power, 2.0, 2.0, 2.0);
  746.  
  747.   return(New);
  748. }
  749.  
  750.  
  751.  
  752. /*****************************************************************************
  753. *
  754. * FUNCTION
  755. *
  756. *   Copy_Superellipsoid
  757. *
  758. * INPUT
  759. *
  760. *   Object - Object
  761. *   
  762. * OUTPUT
  763. *   
  764. * RETURNS
  765. *
  766. *   void * - New superellipsoid
  767. *   
  768. * AUTHOR
  769. *
  770. *   Dieter Bayer
  771. *   
  772. * DESCRIPTION
  773. *
  774. *   Copy a superellipsoid structure.
  775. *
  776. * CHANGES
  777. *
  778. *   Oct 1994 : Creation.
  779. *
  780. ******************************************************************************/
  781.  
  782. static void *Copy_Superellipsoid(Object)
  783. OBJECT *Object;
  784. {
  785.   SUPERELLIPSOID *New, *Superellipsoid = (SUPERELLIPSOID *)Object;
  786.  
  787.   New = Create_Superellipsoid();
  788.  
  789.   /* Get rid of the transformation created in Create_Superellipsoid(). */
  790.  
  791.   Destroy_Transform(New->Trans);
  792.  
  793.   /* Copy superellipsoid. */
  794.  
  795.   *New = *Superellipsoid;
  796.  
  797.   New->Trans = Copy_Transform(Superellipsoid->Trans);
  798.  
  799.   return(New);
  800. }
  801.  
  802.  
  803.  
  804. /*****************************************************************************
  805. *
  806. * FUNCTION
  807. *
  808. *   Destroy_Superellipsoid
  809. *
  810. * INPUT
  811. *
  812. *   Object - Object
  813. *   
  814. * OUTPUT
  815. *
  816. *   Object
  817. *   
  818. * RETURNS
  819. *   
  820. * AUTHOR
  821. *
  822. *   Dieter Bayer
  823. *   
  824. * DESCRIPTION
  825. *
  826. *   Destroy a superellipsoid.
  827. *
  828. * CHANGES
  829. *
  830. *   Oct 1994 : Creation.
  831. *
  832. ******************************************************************************/
  833.  
  834. static void Destroy_Superellipsoid(Object)
  835. OBJECT *Object;
  836. {
  837.   SUPERELLIPSOID *Superellipsoid = (SUPERELLIPSOID *)Object;
  838.  
  839.   Destroy_Transform(Superellipsoid->Trans);
  840.  
  841.   POV_FREE (Object);
  842. }
  843.  
  844.  
  845.  
  846. /*****************************************************************************
  847. *
  848. * FUNCTION
  849. *
  850. *   Compute_Superellipsoid_BBox
  851. *
  852. * INPUT
  853. *
  854. *   Superellipsoid - Superellipsoid
  855. *   
  856. * OUTPUT
  857. *
  858. *   Superellipsoid
  859. *   
  860. * RETURNS
  861. *   
  862. * AUTHOR
  863. *
  864. *   Dieter Bayer
  865. *   
  866. * DESCRIPTION
  867. *
  868. *   Calculate the bounding box of a superellipsoid.
  869. *
  870. * CHANGES
  871. *
  872. *   Oct 1994 : Creation.
  873. *
  874. ******************************************************************************/
  875.  
  876. void Compute_Superellipsoid_BBox(Superellipsoid)
  877. SUPERELLIPSOID *Superellipsoid;
  878. {
  879.   Make_BBox(Superellipsoid->BBox, -1.0001, -1.0001, -1.0001, 2.0002, 2.0002, 2.0002);
  880.  
  881.   Recompute_BBox(&Superellipsoid->BBox, Superellipsoid->Trans);
  882. }
  883.  
  884.  
  885.  
  886. /*****************************************************************************
  887. *
  888. * FUNCTION
  889. *
  890. *   intersect_box
  891. *
  892. * INPUT
  893. *
  894. *   P, D       - Ray origin and direction
  895. *   dmin, dmax - Intersection depths
  896. *   
  897. * OUTPUT
  898. *
  899. *   dmin, dmax
  900. *   
  901. * RETURNS
  902. *
  903. *   int - TRUE, if hit
  904. *   
  905. * AUTHOR
  906. *
  907. *   Dieter Bayer
  908. *   
  909. * DESCRIPTION
  910. *
  911. *   Intersect a ray with an axis aligned unit box.
  912. *
  913. * CHANGES
  914. *
  915. *   Oct 1994 : Creation.
  916. *
  917. ******************************************************************************/
  918.  
  919. static int intersect_box(P, D, dmin, dmax)
  920. VECTOR P, D;
  921. DBL *dmin, *dmax;
  922. {
  923.   DBL tmin = 0.0, tmax = 0.0;
  924.  
  925.   /* Left/right. */
  926.  
  927.   if (fabs(D[X]) > EPSILON)
  928.   {
  929.     if (D[X] > EPSILON)
  930.     {
  931.       *dmin = (MIN_VALUE - P[X]) / D[X];
  932.  
  933.       *dmax = (MAX_VALUE - P[X]) / D[X];
  934.  
  935.       if (*dmax < EPSILON) return(FALSE);
  936.     }
  937.     else
  938.     {
  939.       *dmax = (MIN_VALUE - P[X]) / D[X];
  940.  
  941.       if (*dmax < EPSILON) return(FALSE);
  942.  
  943.       *dmin = (MAX_VALUE - P[X]) / D[X];
  944.     }
  945.  
  946.     if (*dmin > *dmax) return(FALSE);
  947.   }
  948.   else
  949.   {
  950.     if ((P[X] < MIN_VALUE) || (P[X] > MAX_VALUE))
  951.     {
  952.       return(FALSE);
  953.     }
  954.  
  955.     *dmin = -BOUND_HUGE;
  956.     *dmax =  BOUND_HUGE;
  957.   }
  958.  
  959.   /* Top/bottom. */
  960.  
  961.   if (fabs(D[Y]) > EPSILON)
  962.   {
  963.     if (D[Y] > EPSILON)
  964.     {
  965.       tmin = (MIN_VALUE - P[Y]) / D[Y];
  966.  
  967.       tmax = (MAX_VALUE - P[Y]) / D[Y];
  968.     }
  969.     else
  970.     {
  971.       tmax = (MIN_VALUE - P[Y]) / D[Y];
  972.  
  973.       tmin = (MAX_VALUE - P[Y]) / D[Y];
  974.     }
  975.  
  976.     if (tmax < *dmax)
  977.     {
  978.       if (tmax < EPSILON) return(FALSE);
  979.  
  980.       if (tmin > *dmin)
  981.       {
  982.         if (tmin > tmax) return(FALSE);
  983.  
  984.         *dmin = tmin;
  985.       }
  986.       else
  987.       {
  988.         if (*dmin > tmax) return(FALSE);
  989.       }
  990.  
  991.       *dmax = tmax;
  992.     }
  993.     else
  994.     {
  995.       if (tmin > *dmin)
  996.       {
  997.         if (tmin > *dmax) return(FALSE);
  998.  
  999.         *dmin = tmin;
  1000.       }
  1001.     }
  1002.   }
  1003.   else
  1004.   {
  1005.     if ((P[Y] < MIN_VALUE) || (P[Y] > MAX_VALUE))
  1006.     {
  1007.       return(FALSE);
  1008.     }
  1009.   }
  1010.  
  1011.   /* Front/back. */
  1012.  
  1013.   if (fabs(D[Z]) > EPSILON)
  1014.   {
  1015.     if (D[Z] > EPSILON)
  1016.     {
  1017.       tmin = (MIN_VALUE - P[Z]) / D[Z];
  1018.  
  1019.       tmax = (MAX_VALUE - P[Z]) / D[Z];
  1020.     }
  1021.     else
  1022.     {
  1023.       tmax = (MIN_VALUE - P[Z]) / D[Z];
  1024.  
  1025.       tmin = (MAX_VALUE - P[Z]) / D[Z];
  1026.     }
  1027.  
  1028.     if (tmax < *dmax)
  1029.     {
  1030.       if (tmax < EPSILON) return(FALSE);
  1031.  
  1032.       if (tmin > *dmin)
  1033.       {
  1034.         if (tmin > tmax) return(FALSE);
  1035.  
  1036.         *dmin = tmin;
  1037.       }
  1038.       else
  1039.       {
  1040.         if (*dmin > tmax) return(FALSE);
  1041.       }
  1042.  
  1043.       *dmax = tmax;
  1044.     }
  1045.     else
  1046.     {
  1047.       if (tmin > *dmin)
  1048.       {
  1049.         if (tmin > *dmax) return(FALSE);
  1050.  
  1051.         *dmin = tmin;
  1052.       }
  1053.     }
  1054.   }
  1055.   else
  1056.   {
  1057.     if ((P[Z] < MIN_VALUE) || (P[Z] > MAX_VALUE))
  1058.     {
  1059.       return(FALSE);
  1060.     }
  1061.   }
  1062.  
  1063.   *dmin = tmin;
  1064.   *dmax = tmax;
  1065.  
  1066.   return(TRUE);
  1067. }
  1068.  
  1069.  
  1070.  
  1071. /*****************************************************************************
  1072. *
  1073. * FUNCTION
  1074. *
  1075. *   evaluate_superellipsoid
  1076. *
  1077. * INPUT
  1078. *
  1079. *   P          - Point to evaluate
  1080. *   Superellipsoid - Pointer to superellipsoid
  1081. *   
  1082. * OUTPUT
  1083. *   
  1084. * RETURNS
  1085. *
  1086. *   DBL
  1087. *   
  1088. * AUTHOR
  1089. *
  1090. *   Dieter Bayer
  1091. *   
  1092. * DESCRIPTION
  1093. *
  1094. *   Get superellipsoid value in the given point.
  1095. *
  1096. * CHANGES
  1097. *
  1098. *   Oct 1994 : Creation.
  1099. *
  1100. ******************************************************************************/
  1101.  
  1102. static DBL evaluate_superellipsoid(P, Superellipsoid)
  1103. VECTOR P;
  1104. SUPERELLIPSOID *Superellipsoid;
  1105. {
  1106.   VECTOR V1;
  1107.  
  1108.   V1[X] = power(fabs(P[X]), Superellipsoid->Power[X]);
  1109.   V1[Y] = power(fabs(P[Y]), Superellipsoid->Power[X]);
  1110.   V1[Z] = power(fabs(P[Z]), Superellipsoid->Power[Z]);
  1111.  
  1112.   return(power(V1[X] + V1[Y], Superellipsoid->Power[Y]) + V1[Z] - 1.0);
  1113. }
  1114.  
  1115.  
  1116.  
  1117. /*****************************************************************************
  1118. *
  1119. * FUNCTION
  1120. *
  1121. *   power
  1122. *
  1123. * INPUT
  1124. *
  1125. *   x - Argument
  1126. *   e - Power
  1127. *   
  1128. * OUTPUT
  1129. *   
  1130. * RETURNS
  1131. *
  1132. *   DBL
  1133. *   
  1134. * AUTHOR
  1135. *
  1136. *   Dieter Bayer
  1137. *   
  1138. * DESCRIPTION
  1139. *
  1140. *   Raise x to the power of e.
  1141. *
  1142. * CHANGES
  1143. *
  1144. *   Oct 1994 : Creation.
  1145. *
  1146. ******************************************************************************/
  1147.  
  1148. static DBL power(x, e)
  1149. DBL x, e;
  1150. {
  1151.   register int i;
  1152.   register DBL b;
  1153.  
  1154.   b = x;
  1155.  
  1156.   i = (int)e;
  1157.  
  1158.   /* Test if we have an integer power. */
  1159.  
  1160.   if (e == (DBL)i)
  1161.   {
  1162.     switch (i)
  1163.     {
  1164.       case 0: return(1.0);
  1165.  
  1166.       case 1: return(b);
  1167.  
  1168.       case 2: return(Sqr(b));
  1169.  
  1170.       case 3: return(Sqr(b) * b);
  1171.  
  1172.       case 4: b *= b; return(Sqr(b));
  1173.  
  1174.       case 5: b *= b; return(Sqr(b) * x);
  1175.  
  1176.       case 6: b *= b; return(Sqr(b) * b);
  1177.  
  1178.       default: return(pow(x, e));
  1179.     }
  1180.   }
  1181.   else
  1182.   {
  1183.     return(pow(x, e));
  1184.   }
  1185. }
  1186.  
  1187.  
  1188.  
  1189. /*****************************************************************************
  1190. *
  1191. * FUNCTION
  1192. *
  1193. *   insert_hit
  1194. *
  1195. * INPUT
  1196. *
  1197. *   Object      - Object
  1198. *   Ray         - Ray
  1199. *   Depth       - Intersection depth
  1200. *   Depth_Stack - Intersection stack
  1201. *   
  1202. * OUTPUT
  1203. *
  1204. *   Depth_Stack
  1205. *   
  1206. * RETURNS
  1207. *
  1208. *   int - TRUE, if intersection is valid
  1209. *   
  1210. * AUTHOR
  1211. *
  1212. *   Dieter Bayer
  1213. *   
  1214. * DESCRIPTION
  1215. *
  1216. *   Push an intersection onto the depth stack if it is valid.
  1217. *
  1218. * CHANGES
  1219. *
  1220. *   Nov 1994 : Creation.
  1221. *
  1222. ******************************************************************************/
  1223.  
  1224. static int insert_hit(Object, Ray, Depth, Depth_Stack)
  1225. OBJECT *Object;
  1226. RAY *Ray;
  1227. DBL Depth;
  1228. ISTACK *Depth_Stack;
  1229. {
  1230.   VECTOR IPoint;
  1231.  
  1232.   if ((Depth > DEPTH_TOLERANCE) && (Depth < Max_Distance))
  1233.   {
  1234.     VEvaluateRay(IPoint, Ray->Initial, Depth, Ray->Direction);
  1235.  
  1236.     if (Point_In_Clip(IPoint, Object->Clip))
  1237.     {
  1238.       push_entry(Depth, IPoint, Object, Depth_Stack);
  1239.  
  1240.       return(TRUE);
  1241.     }
  1242.   }
  1243.  
  1244.   return(FALSE);
  1245. }
  1246.  
  1247.  
  1248.  
  1249. /*****************************************************************************
  1250. *
  1251. * FUNCTION
  1252. *
  1253. *   compdists
  1254. *
  1255. * INPUT
  1256. *   
  1257. * OUTPUT
  1258. *   
  1259. * RETURNS
  1260. *   
  1261. * AUTHOR
  1262. *
  1263. *   Alexander Enzmann
  1264. *   
  1265. * DESCRIPTION
  1266. *
  1267. *   Compare two slabs.
  1268. *
  1269. * CHANGES
  1270. *
  1271. *   Nov 1994 : Creation.
  1272. *
  1273. ******************************************************************************/
  1274.  
  1275. static int compdists(in_a, in_b)
  1276.  CONST void *in_a;
  1277.  CONST void *in_b;
  1278. {
  1279.   DBL a, b;
  1280.  
  1281.   a = *((DBL *)in_a);
  1282.   b = *((DBL *)in_b);
  1283.  
  1284.   if (a < b)
  1285.   {
  1286.     return(-1);
  1287.   }
  1288.  
  1289.   if (a == b)
  1290.   {
  1291.     return(0);
  1292.   }
  1293.   else
  1294.   {
  1295.     return(1);
  1296.   }
  1297. }
  1298.  
  1299.  
  1300.  
  1301. /*****************************************************************************
  1302. *
  1303. * FUNCTION
  1304. *
  1305. *   find_ray_plane_points
  1306. *
  1307. * INPUT
  1308. *   
  1309. * OUTPUT
  1310. *   
  1311. * RETURNS
  1312. *   
  1313. * AUTHOR
  1314. *
  1315. *   Alexander Enzmann
  1316. *   
  1317. * DESCRIPTION
  1318. *
  1319. *   Find all the places where the ray intersects the set of
  1320. *   subdividing planes through the superquadric.  Return the
  1321. *   number of valid hits (within the bounding box).
  1322. *
  1323. * CHANGES
  1324. *
  1325. *   Nov 1994 : Creation.
  1326. *
  1327. ******************************************************************************/
  1328.  
  1329. static int find_ray_plane_points(P, D, cnt, dists, mindist, maxdist)
  1330. VECTOR P, D;
  1331. int cnt;
  1332. DBL *dists, mindist, maxdist;
  1333. {
  1334.   int i;
  1335.   DBL t, d;
  1336.  
  1337.   /* Since min and max dist are the distance to two of the bounding planes
  1338.      we are considering, there is a high probablity of missing them due to
  1339.      round off error. Therefore we adjust min and max. */
  1340.  
  1341.   t = EPSILON * (maxdist - mindist);
  1342.  
  1343.   mindist -= t;
  1344.   maxdist += t;
  1345.  
  1346.   /* Check the sets of planes that cut apart the superquadric. */
  1347.  
  1348.   for (i = 0; i < PLANECOUNT; i++)
  1349.   {
  1350.     d = (D[0] * planes[i][0] + D[1] * planes[i][1] + D[2] * planes[i][2]);
  1351.  
  1352.     if (fabs(d) < EPSILON)
  1353.     {
  1354.       /* Can't possibly get a hit for this combination of ray and plane. */
  1355.  
  1356.       continue;
  1357.     }
  1358.  
  1359.     t = (planes[i][3] - (P[0] * planes[i][0] + P[1] * planes[i][1] + P[2] * planes[i][2])) / d;
  1360.  
  1361.     if ((t >= mindist) && (t <= maxdist))
  1362.     {
  1363.       dists[cnt++] = t;
  1364.     }
  1365.   }
  1366.  
  1367.   /* Sort the results for further processing. */
  1368.  
  1369.   QSORT((void *)(dists), (size_t)cnt, sizeof(DBL), compdists);
  1370.  
  1371.   return(cnt);
  1372. }
  1373.  
  1374.  
  1375.  
  1376. /*****************************************************************************
  1377. *
  1378. * FUNCTION
  1379. *
  1380. *   solve_hit1
  1381. *
  1382. * INPUT
  1383. *   
  1384. * OUTPUT
  1385. *   
  1386. * RETURNS
  1387. *   
  1388. * AUTHOR
  1389. *
  1390. *   Alexander Enzmann
  1391. *   
  1392. * DESCRIPTION
  1393. *
  1394. *   Home in on the root of a superquadric using a combination of
  1395. *   secant and bisection methods.  This routine requires that
  1396. *   the sign of the function be different at P0 and P1, it will
  1397. *   fail drastically if this isn't the case.
  1398. *
  1399. * CHANGES
  1400. *
  1401. *   Nov 1994 : Creation.
  1402. *
  1403. ******************************************************************************/
  1404.  
  1405. static void solve_hit1(Super, v0, tP0, v1, tP1, P)
  1406. SUPERELLIPSOID *Super;
  1407. DBL v0, v1;
  1408. VECTOR tP0, tP1, P;
  1409. {
  1410.   int i;
  1411.   DBL x, v2, v3;
  1412.   VECTOR P0, P1, P2, P3;
  1413.  
  1414.   Assign_Vector(P0, tP0);
  1415.   Assign_Vector(P1, tP1);
  1416.  
  1417.   /* The sign of v0 and v1 changes between P0 and P1, this
  1418.      means there is an intersection point in there somewhere. */
  1419.  
  1420.   for (i = 0; i < MAX_ITERATIONS; i++)
  1421.   {
  1422.     if (fabs(v0) < ZERO_TOLERANCE)
  1423.     {
  1424.       /* Near point is close enough to an intersection - just use it. */
  1425.  
  1426.       Assign_Vector(P, P0);
  1427.  
  1428.       break;
  1429.     }
  1430.  
  1431.     if (fabs(v1) < ZERO_TOLERANCE)
  1432.     {
  1433.       /* Far point is close enough to an intersection. */
  1434.  
  1435.       Assign_Vector(P, P1);
  1436.  
  1437.       break;
  1438.     }
  1439.  
  1440.     /* Look at the chord connecting P0 and P1. */
  1441.  
  1442.     /* Assume a line between the points. */
  1443.  
  1444.     x = fabs(v0) / fabs(v1 - v0);
  1445.  
  1446.     VSub(P2, P1, P0);
  1447.  
  1448.     VAddScaled(P2, P0, x, P2);
  1449.  
  1450.     v2 = evaluate_superellipsoid(P2, Super);
  1451.  
  1452.     /* Look at the midpoint between P0 and P1. */
  1453.  
  1454.     VSub(P3, P1, P0);
  1455.  
  1456.     VAddScaled(P3, P0, 0.5, P3);
  1457.  
  1458.     v3 = evaluate_superellipsoid(P3, Super);
  1459.  
  1460.     if (v0 * v2 > 0.0)
  1461.     {
  1462.       if (v1 * v2 > 0.0)
  1463.       {
  1464.         /* This should be impossible, since v0 and v1 were opposite signs,
  1465.            v2 must be either 0 or opposite in sign to either v0 or v1. */
  1466.  
  1467.         Error("internal failure in function solve_sq_hit1: %d, %g, %g, %g", i, v0, v1, v2);
  1468.       }
  1469.       else
  1470.       {
  1471.         if (v0 * v3 > 0.0)
  1472.         {
  1473.           if (x < 0.5)
  1474.           {
  1475.             Assign_Vector(P0, P3);
  1476.           }
  1477.           else
  1478.           {
  1479.             Assign_Vector(P0, P2);
  1480.           }
  1481.         }
  1482.         else
  1483.         {
  1484.           /* We can move both ends. */
  1485.  
  1486.           Assign_Vector(P0, P2);
  1487.  
  1488.           Assign_Vector(P1, P3);
  1489.         }
  1490.       }
  1491.     }
  1492.     else
  1493.     {
  1494.       if (v0 * v3 > 0.0)
  1495.       {
  1496.         /* We can move both ends. */
  1497.  
  1498.         Assign_Vector(P0, P3);
  1499.  
  1500.         Assign_Vector(P1, P2);
  1501.       }
  1502.       else
  1503.       {
  1504.         if (x < 0.5)
  1505.         {
  1506.           Assign_Vector(P1, P2);
  1507.         }
  1508.         else
  1509.         {
  1510.           Assign_Vector(P1, P3);
  1511.         }
  1512.       }
  1513.     }
  1514.   }
  1515.  
  1516.   if (i == MAX_ITERATIONS)
  1517.   {
  1518.     /* The loop never quite closed in on the result - just use the point
  1519.        closest to zero.  This really shouldn't happen since the max number
  1520.        of iterations is enough to converge with straight bisection.  */
  1521.  
  1522.     if (fabs(v0) < fabs(v1))
  1523.     {
  1524.       Assign_Vector(P, P0);
  1525.     }
  1526.     else
  1527.     {
  1528.       Assign_Vector(P, P1);
  1529.     }
  1530.   }
  1531. }
  1532.  
  1533.  
  1534.  
  1535.  
  1536. /*****************************************************************************
  1537. *
  1538. * FUNCTION
  1539. *
  1540. *   check_hit2
  1541. *
  1542. * INPUT
  1543. *   
  1544. * OUTPUT
  1545. *   
  1546. * RETURNS
  1547. *   
  1548. * AUTHOR
  1549. *
  1550. *   Alexander Enzmann
  1551. *   
  1552. * DESCRIPTION
  1553. *
  1554. *   Try to find the root of a superquadric using Newtons method.
  1555. *
  1556. * CHANGES
  1557. *
  1558. *   Nov 1994 : Creation.
  1559. *
  1560. ******************************************************************************/
  1561.  
  1562. static int check_hit2(Super, P, D, t0, P0, v0, t1, t, Q)
  1563. SUPERELLIPSOID *Super;
  1564. DBL t0, t1, v0, *t;
  1565. VECTOR P, D, P0, Q;
  1566. {
  1567.   int i;
  1568.   DBL dt0, dt1, v1, deltat, maxdelta;
  1569.   VECTOR P1;
  1570.  
  1571.   dt0 = t0;
  1572.   dt1 = t0 + 0.0001 * (t1 - t0);
  1573.  
  1574.   maxdelta = t1 - t0;
  1575.  
  1576.   for (i = 0; (dt0 < t1) && (i < MAX_ITERATIONS); i++)
  1577.   {
  1578.     VEvaluateRay(P1, P, dt1, D)
  1579.  
  1580.     v1 = evaluate_superellipsoid(P1, Super);
  1581.  
  1582.     if (v0 * v1 < 0)
  1583.     {
  1584.       /* Found a crossing point, go back and use normal root solving. */
  1585.  
  1586.       solve_hit1(Super, v0, P0, v1, P1, Q);
  1587.  
  1588.       VSub(P0, Q, P);
  1589.  
  1590.       VLength(*t, P0);
  1591.  
  1592.       return(TRUE);
  1593.     }
  1594.     else
  1595.     {
  1596.       if (fabs(v1) < ZERO_TOLERANCE)
  1597.       {
  1598.          VEvaluateRay(Q, P, dt1, D)
  1599.  
  1600.          *t = dt1;
  1601.  
  1602.          return(TRUE);
  1603.       }
  1604.       else
  1605.       {
  1606.         if (((v0 > 0.0) && (v1 > v0)) || ((v0 < 0.0) && (v1 < v0)))
  1607.         {
  1608.           /* We definitely failed */
  1609.  
  1610.           break;
  1611.         }
  1612.         else
  1613.         {
  1614.           if (v1 == v0)
  1615.           {
  1616.             break;
  1617.           }
  1618.           else
  1619.           {
  1620.             deltat = v1 * (dt1 - dt0) / (v1 - v0);
  1621.           }
  1622.         }
  1623.       }
  1624.     }
  1625.  
  1626.     if (fabs(deltat) > maxdelta)
  1627.     {
  1628.        break;
  1629.     }
  1630.  
  1631.     v0 = v1;
  1632.  
  1633.     dt0 = dt1;
  1634.  
  1635.     dt1 -= deltat;
  1636.   }
  1637.  
  1638.   return(FALSE);
  1639. }
  1640.